# Run preamble
source("code/preamble.R")

# Load data
study1_data <- readRDS("data/study1_data.rds")

# Select needed variables
bound_data <- study1_data %>% select(
  rep_mean, rep_sd, dem_mean, dem_sd, study, policy, pid2
)

##################
### PLOT MEANS ###
##################

# Reshaping data
mean_data <- bound_data %>%
  select(pid2, policy, rep_mean, dem_mean, study) %>%
  pivot_longer(
    cols = matches("_mean"),
    names_to = "which_partys_dist",
    values_to = "mean"
  ) %>%
  mutate(which_partys_dist = factor(which_partys_dist,
    levels = c("dem_mean", "rep_mean"),
    labels = c("Democrat", "Republican")
  )) %>%
  mutate(pid2 = case_when(
    pid2 == "Democrat" ~ "Among Democrats",
    pid2 == "Republican" ~ "Among Republicans"
  )) %>%
  group_by(pid2, which_partys_dist, policy, study) %>%
  mutate(`Average Perception` = mean(mean)) %>%
  ungroup()

# Grab actual means
actual_means <- study1_data %>%
  select(pid2, policy, par_pos) %>%
  group_by(pid2, policy) %>%
  summarise(`Reality` = mean(par_pos)) %>%
  rename(which_partys_dist = pid2)

# Merge data
mean_data <- mean_data %>%
  left_join(actual_means, by = c("which_partys_dist", "policy"))

# Reshape data for plot
meanPlot_data <- mean_data %>%
  pivot_longer(
    cols = c(
      `Average Perception`,
      `Reality`
    ),
    names_to = "type",
    values_to = "grand_mean"
  ) %>%
  mutate(type = type %>% factor(levels = c("Reality", "Average Perception")))

# Plot
mean_plot <- meanPlot_data %>%
  ggplot() +
  geom_density(
    aes(
      x = mean, fill = which_partys_dist, group = interaction(which_partys_dist, study)
    ),
    alpha = 0.15,
    color = NA
  ) +
  geom_vline(
    aes(
      xintercept = grand_mean,
      group = interaction(type, which_partys_dist, study),
      linetype = type,
      color = which_partys_dist
    ),
    show.legend = FALSE
  ) +
  scale_fill_manual(values = c("dodgerblue1", "firebrick1")) +
  scale_color_manual(values = c("dodgerblue3", "firebrick3")) +
  labs(x = "Perceived Attitude of Average Partisan (Average of Perceived Distribution)") +
  facet_grid(cols = vars(policy), rows = vars(pid2)) +
  theme_bw() +
  theme(
    axis.title.y = element_blank(),
    strip.background = element_blank(),
    strip.text.x = element_text(size = 10),
    strip.text.y = element_text(size = 10),
    axis.text = element_text(size = 10, color = "black"),
    axis.title.x = element_text(size = 12)
  ) +
  guides(fill = "none")

################
### PLOT SDS ###
################

# Reshaping data
sd_data <- bound_data %>%
  select(pid2, policy, rep_sd, dem_sd, study) %>%
  pivot_longer(
    cols = matches("_sd"),
    names_to = "which_partys_dist",
    values_to = "sd"
  ) %>%
  mutate(which_partys_dist = factor(which_partys_dist,
    levels = c("dem_sd", "rep_sd"),
    labels = c("Democrat", "Republican")
  )) %>%
  mutate(pid2 = case_when(
    pid2 == "Democrat" ~ "Among Democrats",
    pid2 == "Republican" ~ "Among Republicans"
  )) %>%
  group_by(pid2, which_partys_dist, policy, study) %>%
  mutate("Average Perception" = mean(sd)) %>%
  ungroup()

# Grab actual SDs
actual_sds <- study1_data %>%
  select(pid2, policy, par_pos) %>%
  group_by(pid2, policy) %>%
  summarise(Reality = sd(par_pos)) %>%
  rename(which_partys_dist = pid2)

# Merge data
sd_data <- sd_data %>%
  left_join(actual_sds, by = c("which_partys_dist", "policy"))

# Reshape data for plot
sdPlot_data <- sd_data %>%
  pivot_longer(
    cols = c(
      `Average Perception`,
      `Reality`
    ),
    names_to = "type",
    values_to = "grand_mean"
  ) %>%
  mutate(type = type %>% factor(levels = c("Reality", "Average Perception")))

# Plot
scaleFUN <- function(x) sprintf("%.2f", x)

sd_plot <- sdPlot_data %>%
  ggplot() +
  geom_density(
    aes(
      x = sd, fill = which_partys_dist, group = interaction(which_partys_dist, study)
    ),
    alpha = 0.15,
    color = NA
  ) +
  geom_vline(
    aes(
      xintercept = grand_mean,
      group = interaction(type, which_partys_dist, study),
      linetype = type,
      color = which_partys_dist
    ),
    show.legend = FALSE
  ) +
  scale_y_continuous(labels = scaleFUN) +
  scale_fill_manual(values = c("dodgerblue1", "firebrick1")) +
  scale_color_manual(values = c("dodgerblue3", "firebrick3"), guide = "none") +
  labs(x = "Perceived Within-Party Attitude Diversity (SD of Perceived Distribution)") +
  facet_grid(cols = vars(policy), rows = vars(pid2)) +
  theme_bw() +
  theme(
    axis.title.y = element_blank(),
    strip.background = element_blank(),
    strip.text.x = element_text(size = 10),
    strip.text.y = element_text(size = 10),
    axis.text = element_text(size = 10, color = "black"),
    axis.title.x = element_text(size = 12)
  ) +
  guides(fill = "none")

##############
### FIGURE ###
##############

# Arrange figure
ggarrange(sd_plot, mean_plot, nrow = 2, labels = c("A", "B"))

# Save figure
ggsave(
  filename = "figures/moment_distributions.png",
  device = "png",
  width = 9,
  height = 7
)

# Print disclaimer
cat(
  "\n\n",
  " See moment_distributions.png for Figure 4."
)

#############################
### CALCULATIONS FOR TEXT ###
#############################

# Underestimating diversity
cat(
  "\n\n",
  " Partisans underestimate the diversity of viewpoints among other partisans:",
  "\n"
)

sd_data %>%
  group_by(pid2, which_partys_dist) %>%
  summarize(
    mean_error_about_sd = mean(sd - Reality),
    se = sd(sd - Reality) / sqrt(n())
  ) %>%
  print()

# Overestimating extremity
cat(
  "\n\n",
  " Partisans generally don't overestimate the extremity of other partisans' viewpoints. However, Republicans do slightly overestimate the extremity of other Republicans' viewpoints.",
  "\n"
)

mean_data %>%
  mutate(
    error_about_extremity = case_when(
      which_partys_dist == "Democrat" ~ Reality - mean,
      which_partys_dist == "Republican" ~ mean - Reality
    )
  ) %>%
  group_by(pid2, which_partys_dist) %>%
  summarize(
    mean_error_about_extremity = mean(error_about_extremity),
    se = sd(error_about_extremity) / sqrt(n())
  ) %>%
  print()

## T-test
t_test <- mean_data %>%
  mutate(
    error_about_extremity = case_when(
      which_partys_dist == "Democrat" ~ Reality - mean,
      which_partys_dist == "Republican" ~ mean - Reality
    )
  ) %>%
  subset(pid2 == "Among Republicans" & which_partys_dist == "Republican") %>%
  pull(error_about_extremity) %>%
  t.test()

cat(
  "\n\n",
  " T-test for whether Republicans overestimate the extremity of other Republicans' viewpoints:",
  "\n",
  " estimate =", t_test$estimate,
  "\n",
  " p =", t_test$p.value
)
